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Using  density  functional  theory  calculations,  many  researchers  have  predicted  that  various 
tungsten-nitride  compounds  WNx(a:  >  1)  will  be  “ultra-incompressible”  or  “superhard”,  i.e.  as 
hard  as  or  harder  than  diamond.  These  compounds  are  predicted  to  have  large  bulk  and  shear 
moduli,  (>  200  GPa)  and  to  be  elastically  and  vibrationally  stable. 

Compounds  with  such  desirable  properties  must  be  energetically  stable  against  decomposition  into 
other  compounds.  This  stability  can  only  be  found  after  the  determination  of  the  convex  hull  for 
WjNi-j  lines  which  connect  the  lowest  enthalpy  structures  as  a  function  of  composition.  The  phase 
diagram  of  the  W-N  structure  is  uncertain,  both  experimentally  and  computationally.  Complex  van 
der  Waals  forces  play  a  significant  role  in  determining  the  structure  of  solid  N2. 

Here  we  use  high-throughput  calculations  to  map  out  the  convex  hull  and  other  low  energy 
structures  for  the  W-N  system.  We  find  that  the  ground  state  of  the  system  is  the  NbO  structure, 
and  that  the  WN2  structures  found  by  Wang  et  al.  are  also  stable  when  van  der  Waals  forces 
are  neglected.  Other  proposed  structures  are  above  the  convex  hull  of  the  W-N  system.  We  show 
how  the  choice  of  density  functional  influences  the  shape  of  the  curve  and  the  structures  that 
form  the  hull.  In  nitrogen-rich  compounds,  the  choice  of  functional  can  dramatically  change  the 
structural  parameters  and  mechanical  behavior.  Using  any  of  the  functionals,  the  bulk  and  shear 
moduli  of  the  NbO  phase  are  comparable  to  the  WNX  compounds  that  have  been  claimed  to  be 
ultra-incompressible  or  superhard. 

PACS  numbers:  64.30.Ef,  64.70.kd,  61.50.Ah 


INTRODUCTION 

The  search  for  materials  with  hardness  comparable  to 
diamond  has  centered  on  compounds  with  highly  direc¬ 
tional,  short,  and  strong  bonds. [1]  In  addition  to  dia¬ 
mond,  these  include  other  proposed  carbon  structures,  [2, 
3]  cubic  boron  nitride[4],  carbonitrides,[5]  and  transition 
metal  borides.  [6,  7]  Over  the  past  several  years  there  has 
been  considerable  theoretical  interest  in  tungsten  nitride 
systems.  [8-11]  These  computational  studies  predict  that 


some  compounds  WNX  with  x  >  2  have  large  bulk  and 
shear  moduli.  In  one  case[ll]  the  hardness  of  the  mate¬ 
rial  was  estimated  to  be  on  the  same  order  as  boro-  and 
carbo-nitrides. 

The  present  studies  all  address  the  issue  of  stability 
for  various  predicted  WNX  compounds.  This  is  done  by 
showing  they  are  elastically  and/or  vibrationally  stable. 
It  is  also  ensured  the  formation  enthalpy  of  the  compound 
is  lower  than  the  formation  energy  of  its  components,  bcc 
tungsten,  and  aN2.[12]  That  is, 


AH(WxNy)  =  [£(WxNy)  -  x  E( W)  -  {y/ 2)  E(aN2)}/(x  +  y)  .  (1) 


where  A//(WxNy)  is  the  formation  energy  per  atom  of 
a  compound  with  stoichiometry  Wa,Ny,  E(Wx’Ny)  is  the 
energy /formula  unit  of  the  compound,  E{ W)  is  the  equi¬ 
librium  energy  per  atom  of  body-centered  cubic  tung¬ 
sten,  and  E(a N2)  is  the  equilibrium  energy  per  molecule 


of  solid  ccN2. 

Condition  (1)  is  necessary,  but  not  sufficient,  for  struc¬ 
tural  stability.  To  be  truly  stable,  a  structure  must 
lie  on  the  convex  hull  constructed  from  the  plot  of  en- 
thalpy/atom  versus  composition  for  all  possible  phases 
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of  the  system.  As  we  shall  see,  most  of  the  WN^  struc¬ 
tures  referred  to  above  do  not  fulfill  this  criterion. 

The  construction  of  the  convex  hull  for  the  tungsten 
nitride  system  is  non-trivial.  The  experimental  phase 
diagram  is  not  well  defined.  [13]  Various  papers  have  pre¬ 
sented  experimental  results  for  cubic[14-17]  (“/3”,  com¬ 
position  33-50%  Nitrogen),  hexagonal  (“5”,  composition 
33-67%  Nitrogen)  [16],  and  tungsten-carbide  (WC)[18] 
structures.  Many  of  these  are  only  found  in  thin  films. 
There  is  no  clear  preference  for  a  ground  state  structure 
at  any  composition.  As  a  result,  there  is  little  guidance 
for  computational  first-principles  studies  of  W-N.  The 
computations  that  do  exist  start  from  a  wide  variety  of 
structures. 

Using  an  evolutionary  technique,  Wang  et  aZ. ,  [9]  found 
that  the  lowest  energy  structure  of  WN2  is  hexagonal. 
The  compound  has  either  an  hP3  structure,  with  space 
group  P6m2,  or  an  hP6  structure  with  space  group 
PGs/mmc.  At  other  stoichiometries  the  calculations 
are  less  certain.  Previous  papers  have  considered  the 
ground  state  of  WN  to  be  NaCl,[19],  tungsten  carbide 
(WC),[8,  10]  or  NiAs.[20]  More  recently,  Song  et  al.  [10] 
looked  at  WN3  in  the  P3Tc  structure,  while  Aydin  et 
aZ.  [11]  have  suggested  that  WN4  has  the  ReP4  structure. 

None  of  these  calculations  addressed  the  complete 
range  of  stoichiometries  in  the  WN  system.  What  is 
needed  is  a  method  which  quickly  examines  a  wide  variety 
of  possible  structures  over  large  ranges  of  composition. 
One  such  method  is  AFLOW,[21]  a  high-throughput  front 
end  for  electronic  structure  calculations.  [22]  AFLOW  al¬ 
lows  us  to  quickly  examine  a  selected  range  of  struc¬ 
tures  using  high-performance  supercomputers  and  mod¬ 
ern  density  functional  electronic  structure  techniques. 
The  AFLOW  prototypes’  database,  which  was  originally 
used  to  describe  intermetallic  alloys,  [23-27]  can  easily 
be  enlarged.  [28]  We  were  therefore  able  to  include  ionic 
and  covalent  structures  which  seem  chemically  similar  to 
W-N.  These  include  borides,  carbides,  oxides,  and  other 
nitrides. 

In  this  paper  we  examine  the  possible  ground  states 
of  the  W-N  system  over  a  wide  range  of  stoichiometries. 
Using  the  power  of  AFLOW,  we  can  examine  hundreds  of 
possible  Wa-Ny  structures.  We  show  that  the  true  ground 
state  of  W-N  is  related  to  what  is  known  as  /3WN,  and 
show  why  it  exists  over  a  wide  range  of  stoichiometries. 
In  addition,  the  bulk  and  shear  moduli  found  for  /3WN 
are  comparable  to  those  WNj  compounds  predicted  to 
be  ultra-incompressible [9]  or  superhard.[ll] 

The  paper  is  organized  as  follows:  Section  describes 
the  computational  methods  used.  Section  describes  the 
main  calculations  of  this  paper:  the  determination  of  the 
convex  hull  for  WN  as  a  function  of  enthalpy  (1).  The 
majority  of  the  calculations  are  done  using  the  Perdew- 
Burke-Ernzerhof  (PBE)[29]  implementation  of  a  Gener¬ 
alized  Gradient  Approximation  density  functional. 

The  ground  state  aN2  structure  of  Nitrogen  is  influ¬ 


enced  by  van  der  Waals  forces.  We  investigate  these 
non-covalent  forces  using  different  density  functionals  in 
Section  .  In  particular,  the  Local  Density  Approxima¬ 
tion  (LDA)  [30,  31]  calculations  are  presented  in  Section  , 
and  the  so-called  DF2  van  der  Waals  functional  (vdW- 
DF2)[32]  is  discussed  in  Section  . 

We  then  discuss  several  specific  structures:  Section 
we  examine  the  cubic  /3WN  phase,  and  the  competing 
<5WN  phase  in  Section  .[16]  Section  looks  at  the  WN2 
structures  found  by  Wang  et  al.,[ 9]  along  with  other  low 
energy  WN2  structures. 

Section  studies  possible  WN3  structures,  including  the 
P3Tc  structure  proposed  by  Song  et  al., [10].  Section 
looks  at  the  ReP4  structure  of  WN4  predicted  by  Aydin 
and  coworkers. 

We  conclude  with  Section  ,  which  summarizes  the  re¬ 
sults,  and  discuss  the  effects  of  the  choice  of  density  func¬ 
tional  on  the  predictions  for  this  system. 

Finally,  as  a  complete  reconstruction  of  a  the  unit  cell 
is  necessary  for  the  reproduction  of  any  calculation,  we 
list  the  crystallographic  information  for  the  most  impor¬ 
tant  structures  discussed  in  this  paper  in  the  supplemen¬ 
tary  material. 

METHODS 

We  began  our  search  for  the  convex  hull  of  the  W-N 
system  by  using  AFLOW  [21,  22,  33]  to  quickly  and  ef¬ 
ficiently  search  through  a  large  database  of  structures. 
As  the  original  AFLOW  prototypes’  database  was  for  bi¬ 
nary  metallic  alloys,  [23-27]  we  extended  it  to  include 
over  fifty  new  structures.  These  include  nitrides,  oxides, 
borides,  and  carbides.  The  important  structures,  includ¬ 
ing  all  of  those  discussed  in  this  paper,  are  described  in 
the  supplementary  material. 

Electronic  structure  calculations  were  done  using  the 
Vienna  Ab  initio  Simulation  Package  (VASP), [34- 
36]  including  core  state  effects  via  the  VASP 
implementation[37]  of  the  Projector  Augmented- 
Wave  (PAW)  method.  [38]  Since  some  our  results  were 
unexpected,  we  checked  them  against  computations 
performed  with  ELK,  an  all-electron  full-potential  Lin¬ 
earized  Augmented  Plane  Wave  (FP-LAPW)  code.  [39] 

AFLOW’s  default  is  to  use  the  Perdew-Burke- 
Ernzerhof  (PBE)  implementation[29]  of  the  Generalized 
Gradient  Approximation  (GGA)  to  Density  Functional 
Theory  (DFT).[40,  41]  There  are,  however,  substantial 
differences  between  PBE  results  and  those  from  the  Local 
Density  Approximation  (LDA). [30,  31]  This  was  shown 
by  Wang  et  aZ.[9],  and  we  shall  see  further  examples  be¬ 
low.  Accordingly,  we  have  done  some  calculations  using 
the  LDA  implemented  in  VASP,  as  shown  in  Section  . 

Molecular  nitrogen,  aN2,  is  a  van  der  Waals  solid,  [42] 
and  neither  the  LDA  nor  the  PBE-GGA  correctly  pre¬ 
dict  its  lattice  constant.  We  studied  the  effect  of  van  der 
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Waals  forces  by  utilizing  the  VASP  implementation[43, 
44]  of  the  vdW-DF2  functional.  [32]  This  functional  was 
developed  to  handle  the  van  der  Waals  interaction  in  gen¬ 
eral  geometries,  so  we  apply  it  over  the  entire  range  of 
compositions.  These  results  are  discussed  in  Section  . 

We  used  the  VASP-supplied  LDA  and  PBE  PAW  po¬ 
tentials  for  Nitrogen  and  Tungsten  (specifically,  May 
2000  LDA  and  April  2002  PBE  PAW  potentials  for  N, 
and  the  July  1998  LDA  and  September  2000  Tungsten 
W_pv  PAWs)  from  the  VASP  distribution.  All  calcula¬ 
tions  use  a  kinetic  energy  cutoff  of  560  eV,  which  is  40% 
larger  than  the  suggested  cutoff  for  Nitrogen  (400  eV). 
Unless  otherwise  stated,  we  let  the  AFLOW  package  de¬ 
termine  the  k-point  mesh  for  each  structure.  This  is  usu¬ 
ally  set  to  use  approximately  the  same  density  of  k-points 
in  reciprocal  space  for  all  structures. 

For  the  NaCl  structure,  e.g,  a  15  x  15  x  15  Monkhorst- 
Pack  grid[45]  was  originally  used.  For  Wang  et  aV s 
P6m2  WN2  structure  we  used  a  15  x  15  x  9  mesh.  In 
the  initial  VASP  calculations  performed  via  AFLOW  all 
structures  are  considered  to  be  spin-polarized.  However, 
in  every  structure  the  self-consistent  magnetic  moment 
was  negligible,  and  the  final  calculations  were  all  done 
assuming  no  moment. 

The  MedeA  software  system[46]  was  used  to  drive 
VASP  in  calculations  of  the  phonon  spectra  and  some 
other  calculations.  Elastic  constants  were  computed  by 
taking  finite  strains  [47,  48]  and  determining  the  slope 
of  the  corresponding  stress-strain  curve  (as  implemented 
in  the  MedeA  package,  or  using  the  native  VASP  op¬ 
tion).  Phonon  frequencies  were  determined  via  the 
frozen-phonon  approximation  using  the  MedeA  package. 

THE  GROUND  STATE  HULL  AND  LOW-LYING 
STATES  OF  TUNGSTEN  NITRIDE 

Overall  we  considered  350  possible  structures  from  the 
expanded  AFLOW  prototypes’  database.  The  equilib¬ 
rium  configurations  for  the  important  structures  are  de¬ 
scribed  in  the  supplementary  material.  Fig.  1  shows  the 
formation  enthalpy  per  atom  compared  to  o;N2  and  bcc- 
W,  as  defined  in  (1). 

The  other  structures  on  the  convex  hull  in  Fig.  1  are 
the  two  WN2  structures  found  by  Wang  et  a£.[9],  and  the 
NbO  structure  (Fig.  2),  which  will  be  considered  further 
in  Section  .  To  our  knowledge  this  structure  has  never 
been  considered  as  a  candidate  structure  for  WN.  Its  ex¬ 
istence  is  not  a  complete  surprise:  WN  has  the  same 
number  of  valence  electrons  as  NbO,  and  the  Nb  and  W 
atoms  have  similar  valence  shells  (4d45s  for  Nb,  5d46s2 
for  W).  The  electron  gained  by  replacing  Nb  with  W  is 
lost  again  as  we  replace  O  with  N. 

The  NbO  structure  is  a  supercell  of  the  NaCl  struc¬ 
ture,  with  ordered  vacancies  on  both  the  Na  and  Cl  sites. 
Many  of  the  other  low-lying  structures  shown  in  Fig.  1 


can  be  constructed  by  removing  atoms  from  a  supercell 
of  the  parent  structure.  Starting  with  the  NaCl  (Bl),  cu¬ 
bic  ZnS  (B3),  and  NiAs  (B81)  structures,  we  constructed 
supercells  of  the  original  cell  (as  described  below)  and 
removed  atoms  from  both  sites  in  a  systematic  manner. 
These  structures  are  color/shape  coded  in  Fig.  1.  In  all 
three  cases,  some  supercells  with  multiple  vacancies  had 
lower  energy  than  the  parent  structure.  The  ground  state 
NbO  can  be  constructed  by  removing  atoms  from  the  par¬ 
ent  NaCl  compound  (NbO),  while  the  P6m2  WN2  struc¬ 
ture  can  be  derived  from  either  the  NiAs  or  the  tungsten 
carbide  structure. 

Some  structures  which  have  been  predicted  to  be  stable 
actually  have  enthalpies  far  from  the  convex  hull.  In 
particular,  NaCl  is  over  0.3  eV/atom  above  NbO,  and 
the  WC  and  NiAs  structures  are  over  0.2  eV/atom  above 
NbO.  We  will  discuss  the  predicted  P3TC  structure  for 
WN3  and  ReP4  structure  for  WN4  in  later  sections. 


RESULTS  USING  DIFFERENT  DENSITY 
FUNCTIONALS 

As  noted  above,  solid  aN2  is  bound  by  the  van 
der  Waals  interaction  between  molecules.  Furthermore, 
Wang  et  al.  found  a  substantial  change  in  the  enthalpy 
of  their  WN2  systems  when  changing  from  PBE  to  LDA 
functionals.  This  suggests  that  the  choice  of  functional 
may  well  change  our  predictions  of  the  structure  in  the 
W-N  system.  We  tested  this  by  comparing  results  for 
selected  structures  using  three  different  density  function¬ 
als:  the  Perdew-Burke-Ernzerhof  Generalized  Gradient 
Expansion  (PBE), [29]  the  Local  Density  Approximation 
(LDA), [30,  31]  and  the  van  der  Waals  functional  devel¬ 
oped  by  Dion  et  al.  (vdWDF2),[32]  as  implemented  in 
VASP  by  Klirnes  et  al.  [43,  44]  We  began  our  tests  by  con¬ 
sidering  the  known  ground  states  of  the  end-points,  aN2 
and  bcc-W. 

First  consider  nitrogen.  There  is  a  computational 
problem  in  using  VASP  to  compute  the  equilibrium  lat¬ 
tice  constant  for  all  of  the  nitrogen  structures  consid¬ 
ered  here.  The  bulk  modulus  of  the  molecular  N2  crystal 
in  any  phase  is  extremely  small,  about  1  GPa.[49]  As  a 
result,  algorithms  which  stop  searching  for  a  minimum 
when  the  pressure  reaches  some  small  number  will  fail 
here  unless  an  extremely  small  tolerance  is  used.  In  addi¬ 
tion,  while  the  calculation  of  the  total  energy  of  a  system 
is  variational,  that  of  the  pressure  is  not.  As  a  result, 
the  pressure  calculation  will  not  converge  unless  we  use 
a  much  larger  plane- wave  basis  set.  We  eliminated  this 
problem  for  the  N2  phases  by  explicitly  calculating  E(V) 
at  discrete  points,  bounding  the  minimum,  and  using  a 
fourth-order  Birch  fit  [50]  to  determine  the  equilibrium 
lattice  constant. 

Experiments  have  determined  that  aN2  is  cubic  with 
either  space  group  P2i3-T4  (#195),  which  has  no  in- 
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FIG.  1:  (Color  online)  Low-enthalpy  structures  of  Tungsten  Nitride,  calculated  using  PBE.  Each  cross  gives  the  enthalpy 
of  a  structure  from  the  AFLOW  database,  as  described  in  the  supplementary  materials.  The  solid  colored  symbols  denote 
structures  constructed  by  removing  atoms  from  a  base  structure  and  allowing  the  system  to  relax,  as  described  in  the  text  and 
the  supplementary  material.  Red  squares:  NaCl  and  constructed  by  removing  atoms  from  supercells  of  NaCl,  including  possible 
/3-phase  structures.  [14]  Green  diamonds:  the  (5-phase  and  related  structures.  [16]  Upward  pointing  purple  triangles:  Cubic  ZnS 
(B3)  and  related  structures.  Downward  pointing  blue  triangles:  NiAs  (B8i).  Gray  pentagons:  WC  and  related  structures.  The 
black  circles  indicate  structures  which  other  researchers  have  described  as  hard. [9,  11]  Note  that  Wang  et  al.'s  P6m2  WN2 
structure  could  be  grouped  with  the  NiAs  structures  as  well  as  the  WC  structure. 


version  site,  or  Pa3-T®  (#205),  which  contains  an 
inversion.  [12]  In  VASP  the  energy  and  structural  differ¬ 
ences  between  the  two  structures  is  negligible.  In  fact,  the 
energy  difference  between  the  five  N2  structures  shown 
in  Fig.  1  is  less  than  4  meV/atom.  This  being  the  case, 
here  we  present  only  the  results  for  the  highest  symmetry 
structure,  Pa3.  The  results  for  the  equilibrium  geometry 
of  this  phase  for  all  three  functionals  are  given  in  Table  I. 

The  PBE  GGA  overestimates  the  equilibrium  lattice 
constant  of  aN2  by  9%  while  the  LDA  underestimates  it 
by  almost  8%,  as  was  also  found  by  Mailhoit  et  aL[52] 
Both  of  these  errors  are  much  larger  than  found  for  most 
elements.  [53]  If  we  consider  the  van  der  Waals  interaction 
using  the  vdW-DF2  functional,  we  get  an  error  in  the 
lattice  constant  of  5.2%,  which,  while  not  ideal,  is  a  sub¬ 
stantial  improvement  on  both  the  LDA  and  PBE  results. 
All  of  the  calculations  give  essentially  the  same  distance 
for  the  N-N  bond,  suggesting  that  errors  in  the  calcula- 


TABLE  I:  Equilibrium  lattice  constant,  internal  parameter, 
and  equilibrium  bulk  modulus  for  the  Pa3  structure  of  QN2, 
as  determined  using  various  density  functionals  and  compared 
to  experiment.  [12,  49]  The  lattice  is  simple  cubic,  and  the 
nitrogen  atoms  sit  on  the  (8c)  Wyckoff  position,  which  has 
one  internal  parameter,  x.  The  equilibrium  lattice  constant 
and  bulk  modulus  are  determined  from  a  fourth  order  Birch 
fit.  The  quantity  d(N-N)  is  the  length  of  the  nitrogen-nitrogen 
bond  in  N2  molecules. 


Functional 

PBE 

LDA 

vdW 

Exp  [12] 

a  (A) 

6.187 

5.223 

5.511 

5.659 

X 

0.052 

0.061 

0.058 

0.056 

d(N-N)  (A) 

1.11 

1.10 

1.11 

1.10 

K0  (GPa) 

0.788 

5.70 

4.69 

1.2  [49] 

tions  are  due  to  the  long-range  interaction  between  the 
molecules.  Both  LDA  and  vdW-DF2  substantially  over- 
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FIG.  3:  The  relative  enthalpy  (1)  of  structures  making  up  the 
ground  state  hull,  and  other  selected  structures  for  the  Tung¬ 
sten  Nitride  system.  These  were  calculated  using  the  Perdew- 
Burke-Ernzerhof  (PBE)  functional [29]  with  VASP.  The  state 
labeled  cI36  at  x  =  1/3  is  the  NbO  structure  with  additional 
defects,  as  described  in  Section  . 


FIG.  2:  (Color  Online)  WN  in  the  NbO  structure,  con¬ 

structed  by  taking  equally  spaced  vacancies  from  both  sites 
of  the  NaCl  structure.  The  PBE  functional  predicts  the  equi¬ 
librium  lattice  constant  to  be  4. 131  A. 


fects  the  predicted  shape  and  structural  composition  of 
the  convex  hull. 


TABLE  II:  Equilibrium  lattice  constant  and  bulk  modulus  for 
bcc  tungsten,  as  determined  using  various  density  functionals 
and  compared  to  experiment. [12,  51]  The  equilibrium  lattice 
constant  and  bulk  modulus  are  computed  from  a  fourth  order 
Birch  fit  to  energy  versus  volume  data. 


Functional 

PBE 

LDA 

vdW 

Exp 

a  (A) 

3.190 

3.143 

3.250 

3.165[12] 

K0  (GPa) 

304 

337 

267 

323[51] 

estimate  the  bulk  modulus,  while  PBE  underestimates  it 
by  about  35%.  The  vdW-DF2  prediction  demonstrates 
the  importance  of  van  der  Waals  interactions  in  Nitrogen 
systems. 

Similar  calculations  were  carried  out  for  the  ground 
state  body-centered  cubic  structure  of  tungsten.  The  re¬ 
sults  are  presented  in  Table  II.  In  this  case  the  equi¬ 
librium  lattice  constant  determined  by  VASP  is  almost 
exactly  identical  to  the  result  from  a  Birch  fit  to  energy- 
volume  data.  We  use  the  latter  so  that  we  can  also  com¬ 
pute  the  equilibrium  bulk  modulus  and  compare  it  to 
experiment.  [51]  All  three  functionals  produce  reasonable 
values  for  the  equilibrium  lattice  constant  and  bulk  mod¬ 
ulus,  with  the  LDA  yielding  perhaps  the  best  agreement. 

The  conclusion  of  this  brief  study  is  that  none  of  the 
three  functionals  is  ideal  for  describing  both  tungsten  and 
nitrogen  structures.  All  three  do  give  the  same  ordering 
of  the  low-energy  states  at  both  endpoints.  As  we  shall 
see  below,  this  agreement  does  not  extend  to  W^Ni-a, 
compounds,  and  the  choice  of  functional  significantly  af- 


LDA  results 

The  calculations  leading  to  the  results  in  Table  I  show 
that  the  PBE  functional  improperly  describes  the  ground 
state  structure  of  aN2.  This  calls  into  question  the  pre¬ 
dicted  shape  of  the  ground  state  hull  shown  in  Fig.  1.  We 
examine  the  possible  errors  in  the  PBE  by  redoing  a  sub¬ 
set  of  our  calculations  with  other  density  functionals.  We 
have  selected  a  data-set  of  34  structures,  including  end 
points,  low-lying  states  (WN2,  NbO  and  related  struc¬ 
tures,  etc.),  “interesting”  states  (P3TC,  ReP4),  and  “par¬ 
ent”  structures  (NaCl,  B3-ZnS).  All  these  are  described 
in  the  supplementary  material.  The  resulting  enthalpy 
calculation  is  shown  in  Fig.  3  for  the  PBE  functional, 
and  Fig.  4  for  LDA. 

The  composition  of  the  hull  is  unchanged  when  go¬ 
ing  from  PBE  to  LDA.  Its  shape  changes  visibly,  as  the 
hexagonal  WN2  phases  are  now  competitive  with  the 
NbO  phase.  Within  the  LDA,  all  structures  are  more 
tightly  bound  compared  to  the  endpoints,  by  a  factor  of 
two  or  more.  On  the  tungsten  right  side  of  the  LDA  plot 
two  structures,  8^  (discussed  in  Section  )  and  MoS2,  are 
extremely  close  to  the  tie-line  between  the  NbO  structure 
and  bcc  tungsten.  This  indicates  that  the  LDA  func¬ 
tional  may  well  be  the  best  choice  for  calculations  on  the 
tungsten-rich  side,  since  several  of  8  phases  are  observed 
experimentally.  [16] 
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FIG.  4:  The  relative  enthalpy  (1)  of  structures  making  up  the 
ground  state  hull,  and  other  other  structures  for  the  Tung¬ 
sten  Nitride  system.  These  were  calculated  using  the  Local 
Density  Approximation  (LDA)  with  VASP.  The  state  labeled 
C.I36  at  x  =  1/3  is  the  NbO  structure  with  additional  defects, 
as  described  in  Section  . 

van  der  Waals  results 

At  first  glance,  Fig.  5,  which  shows  the  vdW-DFT  con¬ 
vex  hull  for  W^Ni-j,  looks  similar  to  the  PBE  (Fig.  1) 
and  LDA  (Fig.  4)  diagrams.  On  closer  examination  we 
see  substantial  differences.  The  principle  difference  is 
that  the  vertex  of  the  hull  at  x  =  1/3  is  not  one  of  the 
WN2  structures  found  by  Wang  et  al..  Rather,  it  is  a 
cubic  structure  we  have  labeled  cI36.  In  fact,  the  Wang 
structures  are  above  a  line  drawn  from  aN2  and  NbO, 
and  so  would  not  be  on  the  hull  even  if  the  cI36  struc¬ 
ture  was  not  present.  The  ReP4  structure  is  now  rather 
low  in  energy,  about  0.1  eV /atom  away  from  the  tie  line. 
We  will  examine  these  structures  in  more  detail  in  later 
sections. 

CUBIC  (d)  WN  STRUCTURES 

Hagg[14,  18]  identified  the  /3  phase  in  the  WN  system 
as  cubic,  with  the  approximate  composition  W2N.  The 
tungsten  atoms  form  a  face-centered  cubic  lattice  and  the 
nitrogen  atoms  are  randomly  distributed,  presumably  at 
the  octahedral  sites.  [16]  Chiu  and  Chuang[54]  grew  thin 
films  of  nominal  composition  WN  using  metallo-organic 
chemical  vapor  deposition  (MOCVD).  They  found  a  cu¬ 
bic  structure  of  composition  WNI;  with  0.7  <  x  <  1.8. 
Where  the  composition  is  near  W2N  they  found  a  lattice 
constant  of  4.125A,  rising  to  4.154A  at  stoichiometry  and 
to  4.172A  as  x  approaches  1.8.  They  state  that  the  tung¬ 
sten  atoms  remain  on  FCC  sites.  For  x  =  1  the  nitrogen 
atoms  fill  the  octahedral  sites.  For  x  <  1  there  are  va¬ 
cancies  on  the  octahedral  sites,  and  for  x  >  1  nitrogens 
populate  both  the  octahedral  and  tetrahedral  sites. 


N  N3W  NW  NW3  W 
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FIG.  5:  The  relative  enthalpy  (1)  of  structures  making  up 
the  ground  state  hull  and  selected  other  structures  for  the 
Tungsten  Nitride  system.  These  were  calculated  using  vdW- 
DF2  van  der  Waals  functional  with  VASP.  The  state  labeled 
cI36  at  x  =  1/3  is  the  NbO  structure  with  additional  defects, 
as  described  in  the  Section  . 

Computationally,  we  find  something  different,  as  seen 
in  Figures  3-5.  All  three  functionals  predict  NbO  to  be 
the  ground  state  structure.  NbO  has  space  group  Pm3m 
(#221).  The  tungsten  atoms  are  on  the  (3c)  Wyckoff 
sites,  and  the  nitrogen  atoms  are  on  the  (3d)  sites  (or 
visa  versa) .  It  can  be  described  as  an  eight-site  supercell 
of  the  sodium  chloride  structure  with  ordered  vacancies 
on  both  the  Na  and  Cl  sites,  as  shown  in  Fig.  2.  Within 
the  PBE-GGA,  we  find  the  equilibrium  lattice  constant 
to  be  4.13lA. 

This  is  in  reasonable  agreement  with  the  experimen¬ 
tal  value  of  4.125A  found  by  Chiu  and  Chaung[54]  and 
reported  in  the  Alloy  Phase  Diagrams. [13]  Chiu  and 
Cliaung  describe  the  site  as  having  all  FCC  and  octa¬ 
hedral  sites  filled,  i.e.  the  NaCl  (Bl)  structure.  Instead 
we  find  vacancies  on  both  sites.  The  minimization  of  the 
actual  Bl  structure  gives  a  lattice  constant  of  4.366A, 
in  close  agreement  with  other  calculations.  [8,  19]  This 
suggests  that  what  Chiu  and  Chaung  were  seeing  was 
actually  the  NbO  structure,  or  something  very  close  to 
it. 

Table  III  lists  the  structural  and  elastic  parameters 
of  WN  in  the  NbO  structure  using  all  three  function¬ 
als.  The  structure  is  elastically  stable,  with  all  of 
the  Born  criteria[55]  satisfied.  We  also  estimate  the 
shear  modulus  using  the  average  of  the  Hashin-Shtrikman 
bounds.  [56,  57]  The  bulk  and  shear  moduli  are  compa¬ 
rable  to  those  found  by  Wang  et  al.  for  the  predicted 
hexagonal  WN2  structures,  the  WN3  structure  studied 
by  Song  et  al, [10]  and  the  WN4  structure  predicted  by 
Aydin  et  al.  [11]  We  have  not  tried  to  predict  the  actual 
hardness  of  the  material,  but  the  cubic  NbO  structure  of 
WN  is  certainly  difficult  to  compress. 

We  also  checked  the  vibrational  stability  of  the  NbO 
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TABLE  III:  Equilibrium  lattice  and  elastic  constants  of  WN 
in  the  NbO  structure,  (Space  Group  Pm3m  #221,  Wyckoff 
positions  (3c)  and  (3d)).  These  were  computed  by  VASP 
using  the  appropriate  PAW  potentials  for  each  exchange- 
correlation  functional.  Elastic  constants  (in  GPa)  were  com¬ 
puted  by  finite  strain  calculations. [47,  48]  The  isotropic  shear 
modulus  G  is  the  average  of  the  Hashin-Shtrikman  bounds 
for  a  cubic  system.  [56,  57] 

Functional  a  (A)  Cn  C12  C44  B  G 

LDA  4.078  884  140  177  388  239 

PBE  4.131  754  126  173  335  220 

vdW-DF2  4.208  655  120  154  298  192 
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FIG.  6:  (Color  Online)  Phonon  frequencies  along  high- 

symmetry  lines  for  the  NbO  structure  of  WN,  found  via  the 
MedeA  package.  [46]  No  imaginary  phonon  frequencies  were 
found,  indicating  that  this  structure  is  stable  against  small 
amplitude  vibrations. 


structure  by  computing  the  phonons  along  high  symme¬ 
try  lines  using  the  MedeA  software  package.  [46]  The  re¬ 
sults  for  the  PBE  calculation  are  shown  in  Fig.  6.  We 
found  no  imaginary  frequencies  over  the  entire  Brillouin 
zone,  leading  us  to  the  conclusion  that  the  NbO  struc¬ 
ture  of  WN  is  at  least  metastable.  Since  this  is  the  lowest 
energy  structure  found,  we  believe  that  this  is  the  stable 
ground  state  structure  of  WN  at  this  composition. 

The  electronic  band  structure  and  density  of  states  of 
the  NbO  structure  of  WN  are  given  in  Figures  7  and  8, 
respectively.  The  structure  is  metallic,  with  the  tungsten 
d  bands  dominating  the  region  near  and  above  the  Fermi 
level. 

The  prediction  of  the  hitherto  unseen  NbO  ground 
state  was  unexpected  (albeit  perhaps  unsurprising  in  ret¬ 
rospect),  and  could  possibly  be  an  artifact  of  either  the 
PAW  potentials  or  the  PBE  exchange-correlation  approx¬ 
imation.  Accordingly,  we  compared  the  energy  of  the 
NbO  phase  of  WN  to  several  of  the  usual  suspects  for 
the  WN  composition.  We  used  two  techniques: 
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FIG.  7:  The  electronic  band  structure  of  WN  in  the  NbO 
structure  at  the  PBE  equilibrium  (a  =  4.131  A). 


e  -  eF  (eV) 

FIG.  8:  (Color  Online)  The  electronic  density  of  states  and 
angular-momentum  decomposed  density  of  states  for  WN  in 
the  NbO  structure.  These  were  computed  via  the  tetrahedron 
method  in  VASP.  The  primary  contributions  are  from  the 
N-p  and  W-d  states,  with  the  W-d  states  dominating  the 
conduction  band. 


1.  All  electron  calculations  using  the  Elk  FP-LAPW 
code[39]  within  the  PBE  GGA.  Elk  does  not  easily 
do  structural  relaxation,  so  we  used  the  equilibrium 
structural  parameters  from  the  VASP  PBE  GGA 
calculations  and  the  same  k-point  mesh.  We  used 
the  supplied  muffin-tin  files  for  W  and  N,  and  set 
the  energy  cutoff  using  a  value  of  RGmal  =  9.0. 

2.  LDA[30,  31]  calculations  using  the  LDA  PAW  po¬ 
tentials  supplied  with  VASP.  In  this  case  we  allowed 
the  structures  to  fully  relax.  All  other  parameters, 
including  the  energy  cutoff,  we  kept  the  same  as  in 
the  PBE  GGA  calculations. 

Table  IV  shows  the  results.  The  energy  differences  be¬ 
tween  the  VASP  and  Elk  PBE  calculations  are  all  less 
than  0.02  eV/atom.  The  VASP  LDA  results  show  the 


TABLE  IV:  Energy  differences  for  various  structures  of  WN 
with  50-50  stoichiometry,  computed  using  VASP  with  PBE 
PAW,  the  Elk  FP-LAPW  code,  and  VASP  with  LDA  PAW. 
The  Elk  calculations  used  the  relaxed  structures  from  the 
PBE  PAW  runs,  while  the  VASP  LDA  calculations  were  fully 
relaxed.  In  each  case  the  energy  of  the  NbO  structure  was  set 
to  zero,  and  we  show  energy  differences  between  the  struc¬ 
tures  in  eV/atom.  In  the  B81  (NiAs)  calculation  the  N  atom 
was  placed  on  the  As  (Wyckoff  position  2c)  site. 


Structure 

VASP  PBE 

Elk 

VASP  LDA 

NbO 

0 

0 

0 

B81  (NiAs) 

0.219 

0.191 

0.124 

Bh  (WC) 

0.315 

0.303 

0.233 

B1  (NaCl) 

0.590 

0.591 

0.603 

B3 (ZnS) 

0.723 

0.735 

0.637 

same  ordering  of  the  structures,  although  the  energy  dif¬ 
ferences  are  larger  (as  much  as  0.1  eV/atom).  We  con¬ 
clude  that  the  NbO  state  is  the  lowest  energy  structure  of 
all  the  50-50  stoichiometries  structures  we  have  studied. 

Since  the  NbO  ground  state  WN  can  be  considered  as 
an  NaCl  structure  with  ordered  defects,  it  is  logical  to  ask 
if  any  other  pattern  of  defects  produce  low  energy  struc¬ 
tures.  As  mentioned  at  the  beginning  of  this  section, 
the  literature[13,  14,  18,  54]  suggests  that  /3WN  is  Ni¬ 
trogen  deficient,  with  composition  approximately  WN0.7. 
The  simplest  Nitrogen  deficient  structure  we  can  make 
is  to  take  one  atom  out  of  an  8  atom  supercell  of  the 
NaCl  structure,  or,  equivalently,  add  one  Tungsten  atom 
back  into  the  NbO  structure.  Experimentally,  the  proto¬ 
type  for  this  structure  is  S3U4,[58]  and  so  we  will  refer 
to  it  by  that  designation.  As  shown  in  Figs.  3-5,  the 
enthalpy/atom  of  W4N3  is  lower  in  the  S3U4  structure 
than  the  enthalpy  of  WN  in  the  NaCl  structure,  but  it  is 
still  significantly  larger  than  the  enthalpy  of  WN  in  the 
NbO  structure.  As  a  test,  we  also  studied  Nitrogen-rich 
W3N4  in  the  S3U4  structure.  We  found  that  its  enthalpy 
was  also  lower  than  NaCl  -  and  even  lower  than  W4N3 
in  the  S3U4  structure  -  but  still  well  above  the  enthalpy 
of  the  NbO  structure. 

We  continued  along  this  path  by  constructing  vacancy 
patterns  in  NaCl  supercells  containing  8,  16,  and  32 
atoms.  We  considered  a  variety  of  possible  unit  cells 
where  a  given  site  in  the  supercell  was  “on”  (occupied) 
or  “off”  (vacancy).  We  then  fully  relaxed  the  unit  cell, 
and  computed  the  enthalpy  according  to  (1).  For  an  N 
atom  cell  there  are  2N  —  1  possible  combinations.  For 
the  larger  cells  we  have  therefore  currently  placed  con¬ 
straints  on  the  search  algorithm  such  that  the  composi¬ 
tion  is  WNX  with  x  £  [1/2, 3/2].  All  of  these  structures 
can  be  considered  as  candidates  for  Hagg’s  /3  phase. 

For  the  8  atom  supercell  (a  simple  cubic  supercell  of 
NaCl)  we  looked  at  all  255  combinations.  We  found 
34  unique  structures,  including  NbO  itself,  CsCl,  S3U4, 
Re03,[59]  cubic  perovskite  (with  formula  unit  NW4  or 


FIG.  9:  (Color  Online)  The  N7W6  supercell.  The  atoms 
labeled  N  and  W  form  the  NbO  structure.  The  atom  labeled 
N1  was  added  at  one  of  the  vacancy  sites,  maintaining  cubic 
symmetry.  Note  the  substantial  relaxation  of  the  W  atoms 
away  from  the  N1  atom.  Aside  from  NbO  this  is  the  lowest 
enthalpy  structure  of  all  those  studied. 


WN4),  monatomic  fee  ( Strukturbericht  symbol  Al),  and 
simple  cubic  (A/,),  as  well  as  a  variety  of  lower  symmetry 
structures. 

The  16  atom  supercell  is  face-centered  cubic,  with  all 
primitive  vectors  doubled  from  the  original  NaCl  cell.  We 
looked  at  26  structures  with  composition  Wm+4Nn+4, 
where  m,n  >  4. 

The  32  atom  supercell  of  the  NaCl  structure  is  body- 
centered  cubic.  In  this  case  we  only  looked  at  structures 
closely  related  to  the  NbO  structure.  Thus  we  examined 
Wi3N12,  which  has  one  of  the  vacant  W  sites  in  the  NbO 
structure  occupied,  and  W12N11,  which  has  an  extra  N 
atom  removed  from  the  NbO  supercell.  We  looked  at 
18  structures  of  varying  composition.  Aside  from  the 
near  50-50  stoichiometric  structures,  the  most  interesting 
one  has  the  composition  N12Wg.  This  structure  will  be 
described  in  detail  in  Sec.  . 

The  low  enthalpy  structures  generated  by  this  process 
are  plotted  as  red  squares  in  Fig.  1.  All  are  closely  related 
to  NbO,  either  by  the  addition  or  subtraction  of  an  ad¬ 
ditional  N  or  W  atom  from  the  16  or  32  atom  supercells. 
For  example,  the  lowest  energy  structure  after  NbO  plot¬ 
ted  on  the  graph  has  composition  N7W6.  We  constructed 
this  structure  starting  with  the  NbO  structure  expanded 
onto  the  16  atom  supercell.  We  then  placed  an  additional 
N  atom  back  onto  one  of  the  vacant  sites,  and  allowed  the 
system  to  relax.  The  cell  is  shown  in  Fig.  9.  There  is  ev¬ 
ery  reason  to  suspect  that  with  increasing  supercell  size 
we  can  find  a  large  number  of  structures  with  enthalpy 
close  to  NbO,  with  varying  compositions. 
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Our  calculations  show  that  low-lying  /3  structures  must 
have  vacancies  on  both  the  N  and  W  sites.  Since  there 
are  a  large  number  of  these  structures  very  close  to  the 
NbO  state,  entropic  effects  must  be  large.  The  observed 
“sodium  chloride”  structure,  with  lattice  constant  4.12- 
4.14A,[13,  15]  is  most  likely  formed  with  approximately 
25%  vacancies  on  each  site.  These  vacancies  will  no  doubt 
be  ordered  locally,  but  when  averaged  over  the  entire 
crystal  will  look  like  the  NaCl  structure.  The  lattice  con¬ 
stant  of  this  phase  will  be  closer  to  NbO  (4.13lA)  than 
pure  NaCl  (4.366A).  The  experimentally  observed  lattice 
constant  of  4.12-4.14A  is  consistent  with  this  conjecture. 


THE  HEXAGONAL  SH  STRUCTURES 

Khitrova  and  Pinsker[16]  detailed  six  hexagonal  u6h ” 
and  rhombohedral  “Sr”  phases  in  the  W-N  system,  with 
compositions  ranging  from  NW2  (Sr)  to  N2W  (Sr1)- 
Four  of  these  structures  were  found  to  have  vacancies 
on  one  of  the  tungsten  sites.  Since  we  were  primar¬ 
ily  interested  in  the  nitrogen-rich  structures  close  to  the 
convex  hull  we  did  not  do  exhaustive  searches  for  or¬ 
dered  vacancies  in  these  structures.  We  did,  however, 
find  that  the  5TH  structure,  taken  without  defects  and 
thus  having  composition  N2W3,  had  the  lowest  energy  of 
all  the  S  structures.  Fig.  1  shows  that  this  structure  is 
within  0.1  eV/atom  of  the  tie- line  within  the  PBE,  and 
the  structure  is  barely  above  the  tie-line  within  the  LDA 
(Fig.  4).  It  is  quite  likely  that  searches  for  ordered  defect 
structures  in  supercells  of  Sj^  will  show  lower  energies. 


PREDICTED  WN2  STRUCTURES 

In  both  the  LDA  and  PBE  the  lowest  lying  states  with 
composition  WN2  are  the  hexagonal  structures  found  by 
Wang  et  al.  [9].  Our  structural  parameters  agree  with 
theirs,  as  shown  in  the  early  columns  of  Table  V.  We  also 
found  the  elastic  constants  of  both  structures,  including 
the  average  of  the  Voigt[60]  and  Reuss[61]  bounds  on  the 
isotropic  bulk  and  shear  moduli.  These  also  agree  with 
Ref.  9,  confirming  that  these  structures  are  candidates 
for  superlrard  WN. 

Both  of  these  structures  can  be  described  as  N2  dimers 
separated  by  Tungsten  planes.  We  might  therefore  ex¬ 
pect  the  addition  of  van  der  Waals  forces  between  the 
dimers  to  be  significant.  Structurally,  as  shown  in  the 
last  column  of  Table  V,  there  is  little  difference  between 
the  vdW-DF2  calculations  and  the  LDA  or  PBE.  The 
hexagonal  unit  cell  is  expanded  slightly  in  both  struc¬ 
tures,  as  is  the  N2  dimer  separation,  but  the  changes  are 
not  substantial.  The  energetic  difference,  is,  however, 
quite  large.  In  the  LDA  (Fig.  4)  and  PBE  (Fig.  3)  these 
WN2  structures  have  enthalpies/atonr  comparable  to  the 
NbO  structure  of  WN. 


TABLE  V:  Equilibrium  structural  parameters  for  the  hexago¬ 
nal  NjW  structures  found  by  Wang  et  al.[ 9]  The  low  symme¬ 
try  hP3  structure  is  in  space  group  P&m2-D\h  (#187),  with 
nitrogen  atoms  at  the  (2g)  Wyckoff  positions  (OOz)  and  the 
tungsten  atom  on  the  (Id)  site  (|||).  In  the  high  symmetry 
liP6  structure  (P63 /mmc-DQh,  #194)  the  nitrogen  atoms  are 
on  (4e)  sites  (OOz),  and  the  tungsten  atoms  are  on  (2d)  sites 
(|||).  For  the  LDA  and  PBE  functionals  we  show  our  results 
(O),  and  the  results  of  Wang  et  al.  (W).  Rjv-jv  and  Rat-w 
are  the  length  of  the  respective  bonds.  The  elastic  constants 
are  in  GPa.  The  isotropic  bulk  (B)  and  shear  (G)  moduli  are 
averages  of  the  Voigt[60]  and  Reuss[61]  bounds. 


Functional 


LDA 


PBE 


vdW-DF2 


Low  symmetry  P6m2,  hP3 


0 

W 

O 

W 

O 

a  (A) 

2.887 

2.887 

2.928 

2.933 

3.000 

c  (A) 

3.877 

3.877 

3.916 

3.918 

3.974 

Z 

0.1804  0.1804 

0.1813  0.1814 

0.1824 

Rat-at  (A) 

1.399 

1.399 

1.421 

1.421 

1.450 

RAT-111  (A) 

2.077 

2.077 

2.104 

2.104 

2.143 

Gn 

638 

654 

573 

588 

493 

C33 

1056 

1082 

954 

973 

827 

C44 

241 

260 

222 

232 

191 

C12 

207 

213 

183 

191 

159 

Ci  3 

230 

248 

193 

206 

157 

B 

396 

412 

351 

255 

297 

G 

246 

255 

226 

231 

195 

High  symmetry  P63 /mmc,  hP6 

O 

W 

O 

W 

O 

a  (A) 

2.893 

2.893 

2.934 

2.939 

3.007 

c  (A) 

7.714 

7.714 

7.790 

7.796 

7.891 

Z 

0.0898 

0.0898 

0.0902 

0.0902 

0.0907 

Rat-at  (A) 

1.392 

1.392 

1.406 

1.406 

1.431 

Rat-ii,  (A) 

2.078 

2.078 

2.105 

2.105 

2.144 

Cn 

633 

642 

568 

579 

486 

C33 

1051 

1078 

952 

973 

827 

C44 

249 

262 

228 

233 

197 

C12 

213 

217 

187 

195 

161 

C13 

237 

252 

199 

211 

197 

B 

399 

411 

353 

364 

298 

G 

245 

252 

225 

228 

195 

Adding  van  der  Waals  forces  changes  the  shape  and 
composition  of  the  ground-state  convex  hull.  The  ground 
state  structure  of  the  system  is  now  the  body-centered 
cubic  structure  shown  in  Fig.  10.  Although  it  is  not 
obvious,  this  structure  was  constructed  from  a  32-atom 
body-centered  cubic  supercell  of  NaCl.  This  results  in 
a  structure  with  space  group  #299.  Nitro¬ 

gen  atoms  are  at  the  (2a),  (6b),  and  (24h)  Wyckoff  po¬ 
sitions  of  the  BCC  lattice,  while  tungsten  atoms  are  on 
the  (8c),  (12d),  and  (12e)  sites.  We  remove  the  atoms  on 
the  (2a),  (6b),  (8c),  and  (12e)  sites,  leaving  12  N  atoms 
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FIG.  10:  (Color  online)  The  cI36  structure  predicted  for  No  W. 
The  basic  structure  looks  the  same  using  LDA,  PBE,  and 
vdW-DF2  functionals.  The  exact  dimensions  are  described 
in  Table  VI.  Note  that  this  is  a  body-centered  cubic  lattice, 
so  the  center  of  the  cube  and  the  cube  corners  are  equivalent 
sites. 


TABLE  VI:  The  cI36  structure  of  WN2,  constructed  from  a 
32  atom  body-centered  cubic  supercell  of  the  NaCl  structure. 
The  space  group  is  Jm3m-0®  (#229),  the  nitrogen  atoms 
occupy  the  (24h)  Wyckoff  positions  (0 yy),  and  the  tungsten 
atoms  occupy  the  (12d)  sites  (0||).  The  first  N-W-N  angle 
is  for  atoms  lying  on  the  surface  of  the  cube  cell.  The  second 


one  of  the  N  atoms  in  the  interior  of  the  cube 

Functional 

LDA  PBE  vdW 

a  (A) 

y 

W-N  bond  (A) 

10.273  10.383  10.490 

0.1454  0.1455  0.1455 

1.840  1.860  1.879 

/ 

Lngles 

W-N-W 
N-W-N  (face) 
N-W-N  (interior) 

108.7°  108.6°  108.5  ° 

108.7°  108.6°  108.5° 

109.9°  109.9°  110.0° 

at  the  (0,  y,  y)(2Ah)  Wyckoff  positions,  and  6  W  atoms 
on  Wyckoff  sites  (01/41/2)  (6 d).  The  nitrogen  atoms  are 
then  allowed  to  relax  away  from  the  y  =  1/4  value  of  the 
original  supercell.  We  have  relaxed  this  structure,  which 
we  will  designate  by  its  Pearson  symbol,  cI36,  using  all 
three  functionals.  The  results  are  shown  in  Table  VI. 

Within  the  LDA  the  cI36  structure  is  not  particularly 
noteworthy.  Its  enthalpy  is  above  the  tie-line,  well  above 
that  of  our  other  WN2  structures,  including  Mo2N,  an¬ 
other  structure  which  can  be  formed  by  removing  atoms 
from  the  NaCl  structure.  With  the  PBE  functional  the 
cI36  structure  is  less  than  0.1  eV/atom  above  the  ground 


FIG.  11:  The  electronic  band  structure  for  the  body-centered 
cubic  cI36  structure  of  WN2  described  in  the  text.  Calcula¬ 
tions  were  done  using  the  vdW-DF2  functional,  including  van 
der  Waals  forces,  within  VASP.  The  top  valence  band  has  a 
very  small  dispersion,  dropping  to  -0.05  eV  below  the  Fermi 
level  at  P.  The  gap  is  0.75  eV,  and  is  direct  at  T.  The  LDA 
and  PBE  calculations  produce  similar  band  structures. 

state  hull,  below  all  of  our  test  structures  save  for  the 
two  from  Wang  et  al. 

Finally,  when  we  use  the  vdW-DF2  functional  the  cI36 
structure  becomes  the  lowest  energy  structure.  This  is 
not  necessarily  because  the  van  der  Waals  interaction 
promotes  this  structure.  Rather,  as  is  seen  in  Fig.  5,  the 
energy  of  the  hexagonal  WN2  structures  moves  up  com¬ 
pared  to  the  aN2  and  bcc  W  endpoints,  even  though  the 
lattice  constants,  atomic  positions,  and  bond  lengths  are 
little  changed.  Since  there  is  no  experimental  evidence 
for  any  of  the  WN2  structures,  and  all  three  function¬ 
als  are  approximations  with  well-known  limitations,  we 
cannot  determine  which  picture  is  correct. 

The  cI36  structure  is  an  insulator.  The  electronic  band 
structure,  obtained  using  the  vdW-DF2  functional,  is 
shown  in  Fig.  11.  The  phase  is  also  elastically  stable,  as  is 
shown  in  Table  VII.  We  computed  the  elastic  constants 
using  the  finite-strain  method,  [47]  fixing  the  primitive 
cell  in  its  strained  position  while  allowing  the  atoms  to 
relax  within  the  cell,  as  allowed  by  symmetry.  Whereas 
cubic  WN,  in  the  NbO  structure,  is  quite  hard,  the  shear 
moduli  of  the  cI36  structure  is  predicted  to  be  very  small 
by  all  functionals. 

Further  work  on  possible  WN2  phases  was  done  by 
Song  and  Wang, [10]  who  studied  WN2  in  the  CoSb2 
structure.  [62]  As  shown  in  Figs.  3-5,  this  structure  is  0.1- 
0.2  eV  above  the  hexagonal  WN2  structures  for  all  three 
functionals.  It  is  therefore  never  a  candidate  for  the  W-N 
ground  state.  We  also  found  two  other  structures  with 
energies  close  to  the  CoSb2  structure:  Mo2N,[63]  which 
is  another  structure  constructed  by  removing  atoms  from 
a  supercell  of  NaCl,  and  AuSn2.[64]  Neither  structure  is 
a  candidate  for  the  ground  state  of  WN2.  Parameters  of 
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TABLE  VII:  Elastic  constants  for  the  cI36  structure  of  WN2, 
computed  using  the  LDA,  PBE,  and  vdW-DF2  function¬ 
als.  Elastic  constants  were  computed  by  the  finite  strain 
method,  [47]  allowing  the  atoms  to  relax  at  each  strain  while 
fixing  the  unit  cell.  Equilibrium  structural  parameters  are 
taken  from  Table  VI.  B  is  the  equilibrium  bulk  modulus,  and 
all  elastic  constants  are  in  GPa.  The  shear  modulus  G  is  an 
average  of  the  Hashin-Shtrikman  bounds. [56,  57] 


Functional 

G11 

C12 

c44 

B 

G11  -  Ci2 

G 

LDA 

119 

77 

13 

91 

42 

15 

PBE 

112 

72 

13 

85 

41 

16 

vdW-DF2 

104 

65 

12 

78 

39 

14 

the  minimum  energy  structures  for  all  of  these  phases  are 
given  in  the  supplementary  material. 

WN3  STRUCTURES 

Song  and  Wang[10]  considered  WN3  in  the  P3Tc 
structure, [65]  and  found  that  it  has  a  negative  enthalpy 
relative  to  aN2  and  bcc-W.  They  showed  that  the  elas¬ 
tic  constants  met  the  Born  criteria  for  long-wavelength 
stability,  [55]  and  so  claim  that  it  is  a  possible  state  of 
WN3.  As  seen  in  Fig.  4,  within  the  LDA  the  P3Tc  does 
have  negative  enthalpy  compared  to  the  end  points.  How¬ 
ever,  it  is  over  0.2  eV/atom  above  the  tie-line  using  the 
LDA,  PBE,  and  vdW-DF2  functionals.  This  means  that 
the  PTc3  structure  cannot  be  an  equilibrium  state  in  the 
W-N  system,  although  it  is  possible  it  may  exist  as  a 
metastable  state. 

We  also  looked  at  a  competing  structure,  Molybdite 
(Mo03).[66]  This  structure  has  the  same  space  group 
( Pnma )  as  P3Tc  and  the  same  occupied  Wyckoff  po¬ 
sitions.  Yet  when  starting  from  different  initial  condi¬ 
tions  it  relaxes  to  a  much  different  unit  cell,  even  though 
the  energy  is  comparable  to  the  PTc3  structure.  The 
minimum-energy  structures  for  both  P3Tc  and  Mo03  are 
given  in  the  supplementary  material. 

wn4  structure 

The  WN4  structure  proposed  by  Aydin  et  al.  [11]  is  the 
most  interesting  candidate  structure  we  studied.  Using 
ReP4  as  its  prototype,  [67]  within  the  LDA  they  predicted 
it  to  have  relatively  large  bulk  and  shear  moduli  (338 
and  198  GPa,  respectively).  In  addition  its  hardness  was 
estimated  to  be  near  that  of  some  of  the  superhard  boro- 
carbides  and  carbonitrides. 

We  investigated  this  structure  in  some  detail.  We  re¬ 
laxed  the  structure  provided  by  the  authors  of  Ref.  11, 
using  the  LDA,  PBE,  and  vdW-DF2  structures.  Within 
the  LDA  we  quickly  found  an  equilibrium  position  close 
to  the  published  result.  When  we  next  computed  the 


FIG.  12:  (Color  online)  Energy  versus  volume  calculations 
for  WN4  in  the  ReP4  structure,  using  VASP  and  the  LDA 
exchange-correlation  functional.  The  original  calculation  in 
the  upper  left-hand  corner  found  a  structure  close  to  that  re¬ 
ported  previously.  [11]  When  the  structure  was  expanded  past 
a  volume  of  420  A3,  the  energy  started  to  drop  rapidly.  Re¬ 
versing  the  expansion  caused  another  shift  in  the  structure, 
where  the  energy  finally  reached  the  values  shown  near  the 
bottom  of  the  graph.  The  labels  (a)  and  (b)  correspond  to 
the  structure  pictures  in  Fig.  13.  All  the  calculations  use  the 
same  space  group,  Pbca-D\\ (#61),  and  the  N  and  W  atoms 
remained  on  (8c)  Wyckoff  positions. 


energy-versus  volume  curve,  in  order  to  determine  the 
bulk  modulus,  we  inadvertently  expanded  the  structure 
past  the  unit  cell  volume  of  420  A3.  The  structure  quickly 
dropped  in  energy,  as  shown  in  Fig.  12.  Compressing  the 
lattice  then  caused  a  transformation  to  yet  another  con¬ 
figuration,  labeled  (b)  in  Fig.  12.  Note  that  at  no  time 
did  the  structure  leave  the  Pbca  space  group  of  ReP4, 
and  the  atoms  remained  on  (8c)  Wyckoff  sites.  Most 
likely  the  finite  jumps  in  volume  we  applied  allowed  for 
large  relaxations  in  the  atomic  positions,  causing  a  jump 
from  one  local  minimum  to  another.  It  is  possible  that 
further  manipulation  of  the  structure  will  result  in  even 
lower  energies.  We  compare  the  original  structure  (a)  to 
the  (current)  minimum  energy  LDA  structure  (b)  in  the 
first  part  of  Fig.  13. 

We  also  relaxed  the  original  structure  of  Aydin  using 
the  PBE  and  vdW-DF2  functionals.  Both  quickly  re¬ 
laxed  away  from  the  original  structure  and  into  the  re¬ 
sults  shown  in  Fig.  13.  Note  that  in  the  vdW-DF2  struc¬ 
ture,  the  system  has  separated  into  layers  of  WN2  and 
N2  molecules.  It  is  unlikely  that  such  a  system  could  be 
very  hard. 

The  minimum  energy  lattice  constants  and  elastic  con¬ 
stants  for  the  structures  shown  in  Fig.  13  are  given  in 
Table  VIII.  In  the  interest  of  saving  space,  we  list  the 
atomic  positions  in  the  supplementary  material.  While 
the  structure  labeled  (a)  has  rather  large  elastic  constants 
and  so  might  be  considered  a  hard  material,  the  other 
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FIG.  13:  Possible  metastable  primitive  cells  of  N4W  in  the  ReP4  structure.  The  distances  are  the  values  of  the  orthorhombic 
lattice  constants  (a,  b,  c)  which  locally  minimize  the  total  energy.  The  LDA  functional  was  used  for  (a)  and  (b),  which  correspond 
the  the  structures  referred  to  in  Fig.  12.  The  minimum  energy  structure  for  the  PBE  functional  is  shown  in  (c),  and  the  minimum 
energy  structure  for  the  vdW-DF2  functional  is  shown  in  (d).  Minimum  energy  lattice  and  Wyckoff  position  parameters  for  all 
structures  are  given  in  the  supplementary  material. 


structures  do  not.  Since  the  (b)  structure  is  lower  in  en¬ 
ergy  than  the  (a)  structure  within  the  LDA,  we  conclude 
that  the  ReP4  structure  of  N4W  is  not  a  hard  material, 
assuming  it  can  be  made.  We  also  note  that  we  can¬ 
not  guarantee  that  any  of  the  structures  studied  here  is 
the  true  minimum  energy  structure  of  ReP4.  There  are 
too  many  possible  local  minima  to  search  through  at  this 
time. 

The  negative  elastic  constants  for  C44  in  the  LDA  (b) 


structure  and  C55  in  the  PBE  and  vdW-DF2  structures 
violate  the  Born  criteria  for  elastic  stability.  We  are  per¬ 
forming  further  calculations  on  these  systems  in  the  hope 
of  finding  still  lower  energy  structures.  Since  we  have  al¬ 
ready  shown  that  the  ReP4  state  is  not  a  viable  candidate 
for  either  a  hard  or  a  stable  phase  in  the  WN  system,  we 
will  defer  discussion  of  these  results. 
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TABLE  VIII:  Lattice  constants  and  elastic  constants  for  the 
WN4  structures  shown  in  Fig.  13.  Elastic  constants  were  com¬ 
puted  within  VASP  using  finite  distortions  of  the  unit  cell, 
and  include  the  internal  relaxation  of  the  atoms  to  the  strain. 
Structure  labels  refer  to  the  structures  shown  in  Fig.  13.  Lat¬ 
tice  constants  are  in  Angstroms,  while  elastic  constants  are  in 
GPa. 


Structure 

(a) 

(b) 

(c) 

(d) 

Functional 

LDA 

LDA 

PBE 

vdW-DF2 

a 

9.093 

9.662 

13.273 

14.116 

b 

5.314  6.279 

4.927 

6.415 

c 

7.330 

7.532 

8.402 

9.542 

C11 

742.1 

86.8 

12.1 

14.6 

C22 

159.2 

96.6 

138.9 

131.3 

c33 

660.1 

232.6 

279.7 

132.8 

Cl4 

23.7 

-11.9 

55.5 

20.5 

C55 

383.0 

7.5 

-3.2 

-4.7 

Ce6 

74.4 

87.8 

2.8 

1.0 

C12 

102.5 

42.3 

4.6 

5.5 

C\  3 

314.6 

58.8 

8.0 

2.9 

C23 

47.3 

82.4 

69.4 

92.0 

DISCUSSION 

Modern  first-principles  electronic  structure  methods 
make  it  relatively  simple  to  determine  the  minimum  en¬ 
ergy  configuration,  band  structure,  and  total  energy  for 
any  reasonably  sized  structure.  Once  that  is  known,  it 
is  straightforward  (although  frequently  tedious)  to  deter¬ 
mine  if  the  structure  is  stable  against  strains  and  vibra¬ 
tions,  and  thus  a  candidate  for  a  stable  or  metastable 
structure. 

It  is  considerably  harder  to  determine  where  a  struc¬ 
ture  fits  energetically  compared  to  other  structures  of 
similar  composition.  In  particular,  we  cannot  determine 
if  a  given  structure  is  stable  with  respect  to  phase  sep¬ 
aration  or  a  structural  phase  transition  until  we  deter¬ 
mine  the  shape  of  the  convex  hull  on  the  enthalpy  versus 
composition  diagram.  High-throughput  methods  such  as 
AFLOW  make  it  simple  to  determine  this  hull.  This  ap¬ 
proach  can  expediently  flag  low-energy  structures  which 
might  be  candidates  for  metastable  states,  as  well  as 
states  which  might  be  found  via  a  pressure  or  temper¬ 
ature  driven  phase  transition. 

This  paper  presents  a  high-throughput  study  of  the 
tungsten-nitrogen  system.  Although  there  have  been 
many  first-principles  calculations  for  this  system,  [8-11, 
19,  20]  there  has  never  before  been  a  study  of  possi¬ 
ble  structures  with  multiple  compositions.  The  high- 
throughput  calculations  were  done  using  AFLOW,  which 
we  enhanced  by  adding  known  nitride,  oxide,  and  boride 
structures  which  had  not  previously  been  included  in 
the  database.  These  calculations  confirmed  the  exper¬ 
imental  observations  that  tungsten  nitride  forms  a  cubic 


system  based  on  the  rock  salt  structure  with  vacancies 
(/3WN).  [14-17] 

In  fact,  we  found  that  the  NbO  structure  (Fig.  2)  is  the 
best  candidate  for  the  ground  state  structure,  while  other 
states  with  varying  patterns  of  vacancies  produce  struc¬ 
tures  nearly  degenerate  with  NbO-WN.  Most  attempts  to 
form  WN  compounds  used  high  temperatures.  It  is  thus 
likely  that  the  experimental  structure  of  WN  will  be  the 
NaCl  structure  with  approximately  25%  random  vacan¬ 
cies  on  both  the  tungsten  and  nitrogen  sites,  with  some 
local  ordering  of  the  vacancies.  This  is  not  quite  in  agree¬ 
ment  with  previous  experiments,  where  it  appears  that 
the  (5WN  has  vacancies  only  on  the  nitrogen  site.  [14-17] 
The  experimental  literature  on  the  WN  system  is  sparse, 
and  more  data  is  needed  to  confirm. 

Establishing  the  WN  convex  hull  allows  us  to  eval¬ 
uate  the  likelihood  of  finding  other  predicted  forms  of 
WiNi-j.  We  find  that  the  P3TC  structure  predicted  by 
Song  and  Wang,  [10]  as  well  as  the  ReP4  structure  pre¬ 
dicted  by  Ay  din  et  al. ,  are  both  well  above  the  convex 
hull  and  thus  unlikely  to  form,  at  least  without  special 
non-equilibrium  processing.  We  did  find  that  the  hexag¬ 
onal  WN2  structures  predicted  by  Wang  et  al.[ 9]  are  sta¬ 
ble,  and  of  comparable  enthalpy  to  the  NbO  phase  -  at 
least  for  calculations  using  the  generalized-gradient  PBE 
functional  and  the  LDA  functional. 

The  problem  of  determining  the  convex  hull  in  WN  is 
exacerbated  by  the  fact  that  the  ground  state  and  other 
low-energy  structures  of  N2  are  van  der  Waals  solids.  As 
shown  in  Table  I,  neither  the  LDA  nor  the  PBE  func¬ 
tionals  adequately  describe  the  ground  state  aN2.  Both 
correctly  describe  the  bond  length  of  the  N-N  dimer, 
but  LDA  substantially  underestimates  the  equilibrium 
lattice  constant  while  PBE  drastically  overestimates  it. 
Since  molecular  nitrogen  is  van  der  Waals  bound,  we  also 
looked  at  the  system  using  the  vdW-DF2  functional  pro¬ 
posed  by  Dion  et  al.  [32],  which  is  readily  implementable 
in  VASP.  [43,  44]  This  produces  a  much  better,  although 
still  imperfect,  lattice  constant  for  aN2. 

While  it  is  not  entirely  clear  that  the  vdW-DF2  func¬ 
tional  can  be  used  to  study  dense  bulk  systems,  we  have 
applied  it  to  this  system  across  all  compositions.  As  seen 
in  Fig.  5,  this  produces  little  change  on  the  tungsten- 
rich  side  of  the  phase  diagram,  but  a  dramatic  change 
on  the  nitrogen-rich  side.  The  vdW-DF2  functional  fa¬ 
vors  the  cubic  cI36  structure  (Fig.  10  and  Table  VI  for 
WN2,  making  it  competitive  in  entlralpy/atonr  with  the 
NbO  structure.  The  hexagonal  WN2  structures,  which 
are  the  ground  state  structures  in  the  LDA  and  PBE, 
move  up  significantly.  This  is  consistent  with  the  limited 
amount  of  experimental  data  we  have:  no  bulk  hexagonal 
WN2  phase  has  been  seen.  Meanwhile  the  ground  state 
NaCl-like  structures,  with  vacancies,  have  been  seen  ex¬ 
perimentally. 

Much  work  remains  to  be  done  on  this  system.  We 
have  essentially  ignored  the  hexagonal  S  phases,  [16] 
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which  are  seen  experimentally  in  thin  films.  This  requires 
yet  another  exhausting  search  for  vacancy  patterns  in  su¬ 
percells,  this  time  starting  from  hexagonal  and  rhombo- 
hedral  unit  cells.  Given  the  relatively  low  energy  posi¬ 
tion  of  the  SJH  phase,  especially  in  the  LDA  calculation,  a 
proper  study  of  vacancy  positions  is  likely  to  change  the 
right-hand  side  of  the  convex  hull. 

We  also  need  to  study  the  effects  of  pressure.  Does 
the  application  of  pressure  appreciably  change  the  phase 
diagram,  especially  on  the  right-hand  side?  Given  the 
extremely  open  nature  of  the  cI36  structure,  even  mod¬ 
est  pressure  may  trigger  a  first-order  phase  transition  to 
another  structure.  We  also  need  to  perform  more  calcula¬ 
tions  to  establish  the  validity  of  the  vdW-DF2  functional 
when  used  in  bulk  systems.  We  will  pursue  these  calcu¬ 
lations  in  the  future. 

In  conclusion,  we  have  used  high-throughput  density 
functional  calculations  to  study  the  tungsten-nitrogen 
phase  diagram,  both  with  and  without  van  cler  Waals 
forces.  Our  calculations  agree  with  experiment  in  that 
the  dominant  structure  of  W-N  is  the  NaCl  structure 
with  vacancies,  and  show  that  this  structure  can  persist 
over  a  large  range  of  compositions.  It  is  unlikely  that 
many  of  the  previously  predicted  nitrogen-rich  phases 
will  ever  be  seen  experimentally,  as  they  are  well  above 
the  ground  state  hull.  The  only  exceptions  may  be  the 
hexagonal  WN2  phases  predicted  by  Wang  et  aL[9]  These 
are  on  the  ground  state  hull  when  we  use  the  LDA  and 
PBE  functionals,  but  not  when  we  include  van  der  Waals 
forces.  In  addition,  the  bulk  and  shear  moduli  of  the  pre¬ 
dicted  cubic  ground  state  has  rather  large  bulk  and  shear 
moduli,  comparable  with  any  of  the  previously  predicted 
ultra-incompressible  structures.  [9-11] 
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